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Cells of the mature afiT cell repertoire arise from the development in the thymus of bone 
marrow precursors (thymocytes). a/JT cell maturation is characterized by the expression of 
thousands of copies of identical af} T cell receptors and the CD4 and/or CD8 co-receptors 
on the surface of thymocytes. The maturation stages of a thymocyte are: (1 ) double nega- 
tive (DN) (TCR-, CD4- and CD8"), (2) double positive (DP) (TCR+, CD4+ and CD8+), and 
(3) single positive (SP) (TCR+, CD4+ or CD8+). Thymic antigen presenting cells provide the 
appropriate micro-architecture for the maturation of thymocytes, which "sense" the signal- 
ing environment via their randomly generated TCRs. Thymic development is characterized 
by (i) an extremely low success rate, and (ii) the selection of a functional and self-tolerantT 
cell repertoire. In this paper, we combine recent experimental data and mathematical mod- 
eling to study the selection events that take place in the thymus after the DN stage. The 
stable steady state of the model for the pre-DP, post-DP, and SP populations is identified 
with the experimentally measured cell counts from 5.5- to 17-week-old mice. We make use 
of residence times in the cortex and the medulla for the different populations, as well as 
recently reported asymmetric death rates for CD4 and CD8 SP thymocytes. We estimate 
that 65.8% of pre-DP thymocytes undergo death by neglect. In the post-DP compartment, 
91.7% undergo death by negative selection, 4.7% become CD4 SP, and 3.6% become 
CD8 SP Death by negative selection in the medulla removes 8.6% of CD4 SP and 32.1 % 
of CD8 SP thymocytes. Approximately 46.3% of CD4 SP and 27% of CD8 SP thymocytes 
divide before dying or exiting the thymus. 
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1. INTRODUCTION 

T cells are a major component of the adaptive immune system that 
play a crucial role in protection against a wide variety of pathogens. 
The T cell receptor (TCR) is generated by somatic recombination 
and has a vast potential to recognize foreign organisms. T cells 
do not recognize pathogens directly, but rather through binding 
pathogen fragments displayed by major histocompatibility com- 
plex (MHC) proteins on the surface of antigen presenting cells 
(APCs). Since MHC molecules are highly polymorphic, useful 
T cells must be selected for in each individual of the species. 
These T cells must have lineage specific effector functions that may 
include direct lysis, production of cytokines, and ability to regulate 
immune responses. Furthermore, some T cells have the potential 
to drive dangerous autoimmune responses ( 1 ). For all of these rea- 
sons, the development of a T cell repertoire is a highly specialized 
and tightly regulated process (2, 3). It takes place in a dedicated 
organ, the thymus, where unique properties of the microenvi- 
ronment ensure the production of functional, yet self-tolerant T 
cells (4-6). 

Multi-potent precursors travel from the bone marrow to the 
thymus through the blood. When they enter the thymus, the pre- 
cursors that commit to the T cell lineage [or canonical early T cell 
progenitors (7) ] , after a 2 week period on average, transition from 



the double negative (DN) stage, where they do not express the co- 
receptors CD4 and CD8, to the double positive (DP) stage, where 
they express both co-receptors (6) . At this stage, the majority of the 
cells have made productive TCR gene rearrangements and express 
a fully formed a/3 TCR on the cell surface. DP cells are located in 
the cortex region of the thymus, where they use their TCR to sur- 
vey self-peptides presented by MHCs on cortical thymic epithelial 
cells (cTECs) (8). DPs that recognize self-peptide-MHC complexes 
with low affinity undergo positive selection, whereas those with 
high affinity are deleted by negative selection (3). Those DP that 
fail to recognize self-peptide-MHC will undergo apoptosis in a 
process referred to as death by neglect. The DP cells that are posi- 
tively selected will then transition to the single positive (SP) stage, 
where they express either the CD4 or CD8 co-receptor, depending 
upon their MHC class specificity (9). MHC class specificity also 
dictates gene expression changes that will ultimately determine 
the effector functions of that T cell: generally, cytotoxicity for CD8 
T cells and cytokine production for CD4 T cells. All positively 
selected cells, whether MHC class I or class II specific, up-regulate 
the chemokine receptor CCR7, which facilitates their migration 
to the medulla, where they undergo further selection events. The 
medulla contains medullary epithelial cells (mTECs) that express 
tissue-restricted antigens regulated by the nuclear factor Aire (10). 
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Exposure to tissue-restricted antigens allows for further deletion 
of T cells specific for self-antigens they may encounter in the 
periphery. Finally, those cells that have been positively selected, 
yet have avoided negative selection, will mature and migrate to the 
periphery (11). 

Previous efforts to develop mathematical models of thymic 
selection have been based on deterministic approaches or cellular 
automata simulations. These studies have shown the importance 
of (i) thymic antigen diversity on the size of the selected T cell 
repertoire (12), (ii) death rates for the more differentiated thy- 
mocyte subsets (13), (iii) thymocyte proliferation and residence 
times (14), (iv) epithelial networks for thymocyte development 
and migration (15), (v) thymocyte competition for antigen (16), 
(vi) self-pMHC complexes expressed on dendritic cells (DCs) 
(17), (vii) receptor-ligand binding affinity (18), and (viii) a sharp 
threshold in TCR-ligand binding affinity that defines the bound- 
ary between negative and positive selection (19). Recent work by 
Ribeiro and Perelson (20) supports the need to develop appropri- 
ate mathematical models to interpret T cell receptor excision cir- 
cles (or TREC) data, which are used to quantify thymic export (20). 
Sinclair et al. in Ref. (21) bring together experimental immunol- 
ogy with mathematical modeling to conclude that CD8 precursor 
thymocytes are more susceptible to death than CD4 precursors. 
This asymmetry in the death rates underlies the experimentally 
observed CD4:CD8 T cell ratio in the periphery. 

Previous experimental studies have tried to determine the num- 
ber of cells going through positive and negative selection in the 
thymus. However, reports estimating the relative number of cells 
undergoing negative selection compared to positive selection have 
been widely variable. Some find that very few cells undergo neg- 
ative selection; others find that two times more cells undergo 
negative selection than positive selection (22-25). In this report, 
we develop a deterministic mathematical model of T cell develop- 
ment in the thymus. Some of us recently published a report where 
we used a novel approach (Birn _/_ Nur77 GFP mice) that allowed 
us to calculate the number of cells undergoing positive and nega- 
tive selection (26). Using previously published data on the relative 
life-span of DP and SP cells, we estimated the hourly rate of both 
positive and negative selection (26). In this manuscript, we make 
use of (i) a subset of this experimental data, and (ii) the asymmet- 
ric death rates observed for CD4 and CD8 precursor thymocytes 
(21), to develop two mathematical models that will enable us to 
estimate selection rates in the cortex and the medulla, and provide 
a quantitative measure for the stringency of thymic selection. The 
first model (see Section 2.1) allows the identification of the fol- 
lowing parameters: DN thymocyte influx into the cortex, pre-DP 
and post-DP death rates, and pre-DP and post-DP differentiation 
rates. Under the assumption of asymmetric death rates for the CD4 
and CD8 SP thymocytes (21), we extend the first model to provide 
estimates for the following medullary rates (see Section 2.2): CD4 
and CD8 SP death, proliferation, and maturation rates. 

2. MATERIALS AND METHODS 

2.1. A FIRST MODEL OF THYMIC DEVELOPMENT AFTER THE DN STAGE 

In this section, we introduce a deterministic model of thymocyte 
development after the DN stage. This first model will be required 
to calibrate the parameter values of the second model introduced 



in Section 2.2. In particular, and as described in Section 3.2, the 
first model allows the identification of parameter values for the 
following rates: <p, m, i± 2 , Pi, and cp 2 . 

This mathematical model makes use of a data set obtained 
from the analysis of eight C57BL/6 wild type and Bim deficient 
mice (average age 9 weeks), that express a Nur77 GFP transgene 
to indicate TCR signal strength experience (26). Flow cytomet- 
ric analysis, as described in that study, used standard markers to 
define various stages of T cell development in the thymus. The 
Nur77 GFP reporter and Bim deficiency were novel modifications 
that allowed us to quantify cells that normally would be deleted 
by strong TCR signaling. In the mathematical model, we consider 
the following thymocyte populations: m, the population of pre- 
selection DP thymocytes (double positive), that are TCRjS 10 ™ and 
CD69 low (26), n 2 , the population of post-selection DP thymo- 
cytes, that are TCR(S+ and CD69 hi s h (26), and n 3 , the population 
of mature SP (single positive) thymocytes. 

We assume that DN thymocytes differentiate to become pre- 
selection DP thymocytes with rate (cells/ day) <p. We further assume 
that after the DN stage, thymocyte cell fate is determined by the 
TCR signal, which a given thymocyte has received. Sinclair et al. 
used CFSE labeling to show that there is no proliferation at the 
post-DP stage (see Figure Al of their manuscript) (21). Stritesky 
et al. looked at proliferation in the post-DP pool with BrdU label- 
ing, and found no evidence (26). We have, thus, only included 
proliferation in the SP thymocyte population (21, 26). The three 
populations, n\,n 2 , and «3, are involved in the following selection 
events in the cortex and the medulla (see Figure 1): 

• 0 —*■ Hi — flux of DN thymocytes into compartment rt\, 

• «i — V «2 - differentiation from pre-DP (n\) to post-DP (n 2 ) 
thymocytes induced by TCR signal, 

• ri\ — V 0 - death by neglect of pre-DP thymocytes due to lack 
of (or weak) TCR signal, 

• n 2 — ^ «3 - differentiation from post-DP (n 2 ) to SP (03) 
thymocytes sustained by intermediate TCR signal, 

• n 2 0 - apoptosis of post-DP (n 2 ) thymocytes due to strong 
TCR signal, 

• H3 — % periphery - exit of SP thymocytes (M3) to the periphery 
(thymic maturation), 

• «3 — > 2 «3 - proliferation of SP thymocytes ( «3 ) in the medulla, 
and 

• 713 —U 0 - apoptosis of SP (M3) thymocytes due to strong TCR 
signal. 

The time evolution of the three populations can be described by 
the following set of ordinary differential equations (ODEs), which 
are based on the selection events described above: 

— = 4>-<p l m - mni , 

at 

dn 2 . . 

-3- = <P\tl\ - <p 2 n 2 - fl 2 n 2 , UJ 
at 

dn-} 

— = <p 2 n 2 - <?3"3 - £i>3«3 + ^3«3 • 
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FIGURE 1 | Thymic development as hypothesized in the first 
model. The flux, <p, represents the differentiation of DNs into 
pre-DPs. Pre-DP thymocytes have two fates: further differentiation 
into the post-DP pool (p,) or death by neglect (/t,). Post-DP 



thymocytes have two fates: further differentiation into the SP pool 
(<p z ) or death by apoptosis (/x z ). Finally, SP thymocytes have three 
fates: maturation and exit into the periphery (tp 2 ), death by apoptosis 
(fi 3 ), or proliferation U 3 ). 



We are interested in studying the steady state of these popu- 
lations, as the experimental data correspond to population cell 
numbers in the three stages (pre-DP, post-DP, and SP) for a steady 
state thymus (26). The steady state of the system of equations 
[equation (1)] is given by: 



<Pl + Ml 



<P2 + M2 



<P3 + M3 - A 3 



(2) 



This unique steady state exists if and only if <pi + /X3 — X3 > 0, 
so that we have n| > 0. In order to study the linear stability of the 
steady state, we calculate A, the Jacobian matrix of equation (1), as 
follows: 



A : 



-{(pi + Ml) 

<Pl 
0 



0 

<<P2 +M2) 
<p 2 



0 
0 

-(<P3 + ^3 



(3) 



A3) 



A is also the Jacobian matrix at the steady state n* = 
(n*,n|,n|), as the system of ODEs [equation (1)] is linear. 
The three eigenvalues of A are given by (as the matrix is lower 
triangular): 

Pi = -Opi+Mi), Pi = -((pi + Hi), Pi = -OP3+M3-A3) ■ 

Therefore, the steady state [equation (2)] is stable, if and only 
if, <pi + /X3 — A3 > 0, which is also the condition for its existence. 
We conclude this section with the analytical solution of the system 
of ODEs [equation (1)], given initial conditions, which provides 
the time evolution of the three thymocyte populations: 



ni (t) = «* + «i(0) e-fa+^i)' , 

«2(0 



r£ + 



-(<0i+MiK 



[(<p 2 + M2) - (<P\ +AH)] 
+ n 2 (0) e -fe+w)f > 



n 3 (0 = n* + 



«i(0) <pi 



[{<Pl + M2) - (Vi +Ml)l 

<P2 



(4) 



-(^i+/ii)t 



+ 



[(<?3 + M3 - ^3) - (<Pl + Ml)] 



[(^3 + M3 - A3) - ((p 2 + M2)] 
+ « 3 (0) e -te+f3-^)' , 

where Mi(0), H2(0), «3(0) represent the initial conditions for the 
thymocyte populations. Note that in the late time limit, that is, 
if t—> +00 and ^3 + /X3 — A3 > 0, then «i(f) — ► «*, «2(f) —*■ "2 
and «3(f) — > ri^, as it is the unique stable steady state. 

2.2. A SECOND MODEL OF THYMIC DEVELOPMENT AFTER THE DN 
STAGE: CD4 AND CD8 SP THYMOCYTES 

As described in Section 2 . 1 , the first deterministic model will allow 
us to calibrate some of the parameters of a more comprehensive 
model, which we now introduce. We subdivide the SP thymocyte 
population in two classes: CD4 SP and CD8 SP thymocytes. This 
is an extension of the model introduced in the previous section, 
and is motivated by the fact that experimentally, SP thymocytes 
express either the CD4 or the CD8 co-receptor. We now have four 
different thymocyte populations to consider: n\ , the population of 
pre-selection DP (double positive) thymocytes, ni, the population 
of post-selection DP thymocytes, M4, the population of mature 
CD4 + SP (single positive) thymocytes, and «8, the population of 
mature CD8+ SP (single positive) thymocytes. 

As described in the previous section, we assume that DN thy- 
mocytes differentiate to become pre-selection DP thymocytes with 
rate (cells/day) <p, and that after the DN stage, thymocyte cell 
fate is determined by the TCR signal, which a given thymocyte 
has received. Thus, the four populations, n\, ri2, «4, and ng, with 
« 3 = n 4 + n g , are involved in the following selection events in the 
cortex and the medulla (see Figure 2): 



0 



«i - flux of DN thymocytes into compartment ri\ 



n\ — > «2 - differentiation from pre-DP {n{] to post-DP (nz) 
thymocytes induced by TCR signal, 
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FIGURE 2 | Thymic development as hypothesized in the second model 

We assume there is a flux, tp, that represents the differentiation of DNs into 
pre-DPs. Pre-DP thymocytes have two fates: further differentiation into the 
post-DP pool {tp,) or death by neglect (/i,). Post-DP thymocytes have three 



fates: further differentiation into the SP pool as CD4 SPs (tp a ), or CD8 SPs {tp a ), 
or death by apoptosis (ii 2 ). Finally, CD4 or CD8 SP thymocytes have three 
fates: maturation and exit into the periphery (| 4 ) or (f a ), death by apoptosis 
(mJ or (/n s ), or proliferation (A 4 ) or U a ). 



• «i — > 0 - death by neglect of pre-DP thymocytes due to lack 
of (or weak) TCR signal, 

• «2 — ^ «4 - differentiation from post-DP (bj) to CD4+ SP (B4) 
sustained by intermediate TCR signal, 

• «2 —> "8 - differentiation from post-DP (bj) to CD8 + SP (Bg) 
sustained by intermediate TCR signal, 

• «2 — ■*■ 0 - apoptosis of post-DP («2) thymocytes due to strong 
TCR signal, 

• «4 —> periphery - exit of CD4+ SP thymocytes (B4) to the 
periphery (thymic maturation), 

• «8 — > periphery - exit of CD8 + SP thymocytes (Bg) to the 
periphery (thymic maturation), 

• «4 — %■ 2 «4 - proliferation of CD4+ SP thymocytes (M4) in the 
medulla, 

• «8 — ^ 2 «8 - proliferation of CD8 + SP thymocytes (ttg) in the 
medulla, 

■ ti4 — > 0 - apoptosis of CD4+SP (714) thymocytes due to strong 
TCR signal, and 

• Mg — > 0 - apoptosis of CD8+SP (ng) thymocytes due to strong 
TCR signal. 

We assume that all model parameters are positive, that is, fi\, 
t-t-2, M4, M8> <Pu <P2> <P4, <Ps, ?4> §8, 4>, ^4, ^8 > 0, and note that 
the parameters and thymocyte populations of the first and second 
model are related by the following equations: 



<P2 = <Pi + fs 
(i 4 « 4 + /z 8 



§4 «4 +?8«8= <P3 tl 3 , 
= (1 3 « 3 , X 4 « 4 + k & n & 



«3 



(5) 



The time evolution of the four populations can be described by 
the following set of ODEs: 



~w 

dt 

d n,\ 
At 

djn 
di 



<p - tpini - (iim , 

<Pl n l - (<P4 + n) n 2 - j"-2«2 . 
(P4?l2 — ^4«4 — (14H4 + X4H.4 ■ 
(f S n 2 - £s«8 - M8«8 + A.8«8 ■ 



(6) 



We are interested in studying the steady state of these popu- 
lations, as the experimental data correspond to population cell 
numbers in the four stages (pre-DP, post-DP, CD4 SP, and CD8 
SP) for a steady state thymus (26). The steady state of the system 
of equations [equation (6)] is given by: 



<P\ + Ml 
nl<P4 

?4 + M4 — ^4 



<P4 + (Pi + (12 



Hs, + Ms — ^8 



(7) 



This unique steady state exists if and only if £4 + (14 — X4 > 0 
and §8 + Ms — As > 0, so that we guarantee n\ > 0 and «g > 0. In 
order to study the linear stability of the steady state, we calculate 
B, the Jacobian matrix of equation (6), as follows: 



(<Pi + Mi) 
n 
0 

1 0 



0 

l +tp S + 11!) 
<Pi 



0 
0 

-(f 4 + M4 ■ 
0 



A, 4 ) 



-(£ 8 + /«8 — A. 8 )y 



(8) 



B is also the Jacobian at the steady state n* = (n*, n|, n|, «g), 
as the system of ODEs is linear. The four eigenvalues of B are 
given by: 

Pi = -{<P\ + IM) , P2 = -in + n + M2) > 

B 3 = -(£4 + (L4 - X4) , B 4 = -(§ 8 + M8 - A.g) . 

Therefore, the steady state equation (6) is stable if and only if 
£4 + (14 — X4 > 0 and ^8 + /xg — As > 0, which is also the condi- 
tion for its existence. We conclude this section with the analytical 
solution of the system of ODEs equation (6), given initial condi- 
tions, which provides the time evolution of the four thymocyte 
populations: 



m(i) = n\ + ni (0) e -(«+^i) f , 

«2(0 



-((fi+Mi)t 



[((P4 + (ps + M2) - (n + vi)] 

+ n 2 (0) e -<w+w+M2)t > 
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"4(f) = < + 



«i(0) (pi 



[((Pi +(p& + M2) ~ ((Pi + Ml)] 



(P4+< 



-(tpi+lil)t 



+ 



[(£4 + M4 - A. 4 ) - (pi + Ml)] 

«2(0)g4 

[(f 4 + M4 - ^4) - ((Pi +(P&+ M2)] 

+ « 4 (0) e -«4+M4-^4)t ; 

«l(0) ipj 



-(^4+^8+/U2)f 



-(pi+/il)f 



+ 



« 8 (r) = «* h 

[((p 4 + (p s + M2) - ((Pi + Mi)] 

<P4 + <P& 

X . 

[(& + Ms - X 8 ) - (pi +Ml)] 

«2(0) ^8 

[(?8 + M8 - *8) - (<?>4 + (P& + M2)] 

+ n 8 (0) e-^+w-^ , 

(9) 

where ni(0), «2(0), M4(0), ns(0) representthe initial conditions for 
the thymocyte populations. Note that in the late time limit, that 
is, if t— > +00 and §4 + /ii — A4 > 0 and §8 + /xg — As > 0, then 
m(f) ->■ n*, m(t) -> nf, /14(f) ->■ n|, and n 8 (f) -> nj, as it is 
the unique stable steady state. 

3. RESULTS 

3.1. PARAMETER ESTIMATION FOR THE FIRST MODEL (MEANS) 

In this section, we make use of previously published experimental 
data (26) that provide thymocyte cell counts for the three subsets 
considered in the first model: pre-DPs, post-DPs, and SPs. The 
original experiments have been carried out for two types of mice: 
wild type mice (Bim +/ ~) and Bim deficient mice (Bim _/_ ). In 
this paper, we will only be considering the wild type experimental 
results. The data will allow us to provide experimental estimates 
for the steady state thymocyte cell counts: n*, rJ. Note that, 
in order to estimate rates (with units of inverse time), thymo- 
cyte cell counts are not enough. Thus, we will make use of the 
additional knowledge provided by experimentally determined res- 
idence times for each population, r „ with i = 1 , 2, 3. If we make use 
of the model (see Section 2.1), the residence time in compartment 
i can be expressed as: 



1 

(Pi + M> 



for ie {1,2,3}. 



The experimental data (see Table 1) correspond to the num- 
ber of cells (thymocytes) at steady state (26), in each of the 
thymic compartments considered in the mathematical models (see 
Sections 2.1 and 2.2), and for eight different mice (j= 1,2, . . ., 8). 

We have made use of the following average residence times in 
each compartment (27-29) 



Ti = 60 h = 2.5 days , T2 
r 3 = 96 h = 4 days . 



16 h = 0.67 days , 



In order to derive estimates for the model parameters, we have 
carried out the following steps: 

1. We make use of the experimentally determined mature SP thy- 
mocyte flux from the medulla to the periphery, which has been 
estimated to be 1-4 x 10 6 cells per day (14, 26, 30). This flux 
corresponds to about 1% of thymocytes leaving the thymus 
every day (30). Given this flux, which we denote by </> out , «j, 



and the fact that ■ 



<^3 nt , we can obtain an estimate for 



(P2>. We have chosen cj> out to be 2.5 x 10 cells per day (14, 30). 

2. Given r^, <z>3, and the fact that r 3 = — \ — , we can obtain an 
estimate for ^3. 

3. Given t\, «*, and the fact that n* = <f> t\ , we can obtain an 
estimate for (p. 

4. We also have «| = <p\ ti n* , which, in principle, allows us to 
estimate (p\ . We make use of linear regression techniques to do 

so (31, 32). 

Let us introduce a\ by the following equation, a\ = p-, and 

make use of the experimental data to write: n^' 1 = a n*' 1 + e', 
for i= 1, 2, . . ., 8. Thus, the squared error is given by: 



E( ai ) = (nT 



«r ) 



We minimize E(a\) with respect to «i, that is J— 



0. 



Solving for fli , we obtain: 



01 



E8 *, 
,= 1 "l 



E8 / *,!-. Z 

i=l («1 ) 



Table 1 | Experimental steady state thymocyte cell counts for the wild type pre-DR post-DR CD4 SR and CD8 SP populations. 



Mouse 


n* (pre-DP) (cells) 


n* (post-DP) (cells) 


n| (SP) (cells) 


n* (SP CD4) (cells) 


n| (SP CD8) (cells) 


1 


82.58 x 10 6 


9.30 x 10 6 


18.36 x 10 6 


13.85 x 10 6 


4.51 x 10 6 


2 


142.19 x 10 6 


19.94 x 10 6 


26.20 x 10 6 


18.73 x 10 6 


7.46 x 10 6 


3 


89.00 x 10 6 


5.98 x 10 6 


15.98 x 10 6 


11.88 x 10 6 


4.10 x 10 6 


4 


29.32 x 10 6 


2.09 x 10 6 


5.61 x 10 6 


4.40 x 10 6 


1.21 x 10 6 


5 


29.32 x 10 6 


2.09 x 10 6 


5.61 x 10 6 


4.40 x 10 6 


1.21 x 10 6 


6 


51.26 x 10 6 


5.93 x 10 6 


9.01 x 10 6 


6.85 x 10 6 


2.16 x 10 6 


7 


64.48 x 10 6 


6.81 x 10 6 


11.64 x 10 6 


9.03 x 10 6 


2.61 x 10 6 


8 


218.94 x 10 6 


15.42 x 10 6 


40.20 x 10 6 


29.46 x 10 6 


10.74 x 10 6 


Mean 


88.39 x 10 6 


8.45 x 10 6 


16.57 x 10 6 


12.33 x 10 6 


4.25 x 10 6 


Standard deviation 


60.11 x 10 6 


5.89 x 10 6 


11.05 x 10 6 


7.94 x 10 6 


3.12 x 10 6 



The bold font highlights the mean and the standard deviation from the individual mice data. 
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Given a\, we can then estimate (fi from the equation 

5. Given Ti, cdi, and the fact that Ti = — \ — , we can obtain an 
estimate for fi\. 

6. We are now left with three remaining parameters: <p2, £i2> and 
A3. Given the experimental constraints on r\, T2, and T3, we 
assume that the average time to proliferate, I/A3, cannot be 
<7 days. Therefore, we consider A3 to be constrained in the 
interval [1/7, r^ 1 ], with time measured in days. We sample 
equally spaced values for A3, and for each value, we compute 
<P2 = "3(w+w jg) _ ratio «2 = % is computed using the 
linear regression method described above (see Figure 3). In 
this way, we obtain an estimate for ^3. We note that the p- 
values for a\ and «2 are given by 7.57 x 1CP 3 and 6.85 x 10~ 3 , 
respectively (both smaller than the significance level a = 0.05). 

7. Given To, a 2, and the fact that x-> = — 7 — , we can obtain an 

z, Ti> " I-12+<P2 

estimate for [X2- 

8. From steps 6 and 7 above, we have generated (a table of) val- 
ues for <f2 and 112, given a fiducial value for A3 in the interval 
[1/7, t 3 ~']. The mice considered in the experimental study are 
5.5-17 weeks old, and their thymus is in steady state (26). Thus, 
we expect that the parameter values can only be accepted if the 
corresponding system of ODEs attains steady state by 3 weeks. 
Therefore, we only accept parameter values which provide thy- 
mocyte cell counts at time t= 21 days that are within one 
standard deviation from the experimentally determined values 
(see Table 1 ) . That is, we impose for the given parameter set that 
the mathematically predicted value «,(t=21 days) belongs to 
the interval n* ± cr,-, with i = 1,2, 3, and where n* is the (exper- 
imental) mean number of cells in compartment i, and cr, is the 
(experimental) standard deviation in compartment i, as given 
in Table 1 . 

We obtain the following parameter values: 

<p = 35.350 x 10 6 cells/day , m = 0.263 day -1 , 

fl 3 = 0.099 day -1 , <p x = 0.137 day -1 , <p 3 = 0.151 day -1 , 

and 

/x 2 € [1.295, 1.443] day -1 , <p 2 € [0.050,0.198] day -1 , 
A 3 £ [0.143,0.250] day -1 . 

These parameters imply the following thymic selection rates: 

3. 1. 1. Death rates 

9.7 x 10 5 cells/h die by neglect in compartment 1 (/U-i «*), 

4.8 x 10 5 cells/h die by negative selection in compartment 2 
(fi2 n 2.)' and 6.9 x 10 4 cells/h die by negative selection in 
compartment 3 {[±3 n|). 

3. 1.2. Differentiation rates 

5.0 x 10 5 cells/h are positively selected in compartment 1, that is, 
become post-DP from pre-DP (<p\ «*), 4.4 x 10 4 cells/h are posi- 
tively selected in compartment 2, that is, become SP from post-DP 
(q>2 fij), and 1.0 x 10 5 cells/h leave compartment 3 to go to the 
periphery (cp 3 
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FIGURE 3 | Linear regression plots for the first model. 



3. 1.3. Proliferation rate 

1.3 x 10 5 cells/h proliferate in compartment 3 (A3 n|). 

We have also computed the stringency of thymic selection, 
which we define as given by the ratio: 

^— i = 6.79% . 

4> 

Finally, we have computed the (per cell) probability to die, given 
that the cell is in compartment i (i= 1, 2, 3), as well as the (per 
cell) probability to proliferate in the medulla. We have obtained: 

Ml M2 
pi = — — — = 65.8% , p2 = — — — = 91.7% , 

Ml + <Pl " M2 + <Pl 

M3 ^3 
p 3 = — = 22.9% , <J3 = = 42.2% . 

M3 + <P3 + A.3 M3 + P3 + ^3 

3.2. PARAMETER ESTIMATION FOR THE SECOND MODEL (MEANS) 

In this section, we make use of previously published experimental 
data (26) that provide thymocyte cell counts for the four subsets 
considered: pre-DPs, post-DPs, CD4 SPs, and CD8 SPs. We only 
make use of the experimental data for the wild type mice. The data 
will allow us to provide experimental estimates for the steady state 
thymocyte cell counts: n*, n\, b|, n\. As described in Section 3.1, 
we also need residence times for each population subset, r;, with 
z= 1, 2, 4, 8. If we make use of the model (see Section 2.2), the 
residence time in compartment can be expressed as: 

r; = , for i 6 {1,2} , and 

<Pi + fJ-i 

Ti = l - , for i € {4, 8} . 

?i + Mi 

Recent experimental data provide support for asymmetric 
death rates in the CD4 and CD8 SP compartments (21). The 
estimated death rates for CD4 and CD8 SP thymocytes are 1 

= 0.04day~ 1 and /zg = 0.11 day -1 . We also make use of the 
estimates derived in Section 3.1 for <f>, /U-i, M2> <pi> and q>2- Finally, 



1 Private communication from Ben Seddon and Andy Yates. 
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the average residence times in each compartment, as described in 
Section 3.1, are given by: 

n = 60 h = 2.5 days , t2 = 16 h = 0.67 days , 
r 4 = 96 h = 4 days , t 8 = 96 h = 4 days . 

In order to derive estimates for the model parameters, we have 
carried out the following steps: 



l 



we can obtain an 



M8+?8 



, we 



1. Given T4, /X4, and the fact that T4 
estimate for £ 4. 

2. In the same way, given rg, /J-g, and the fact that tg 
can obtain an estimate for 

3. We are now left with four remaining parameters: q>4, (fg, A4, and 
kg. We know that cpi = <P4 + (ft,- We sample ^4 in the interval 
[0, <p{\, where cp2 is the mean value of the interval obtained in 
Section 2.1, and for each fiducial value for ^4, we compute the 
corresponding value for ^8- 

4. Given r 4 , ^4, and the fact that «| = "\ ^ , we can compute 

the fraction 03 = ^ by linear regression (see Figure 4), and 
thus obtain an estimate for A4. Note that we will reject values 
of A.4 that imply the proliferation time is larger than 7 days (see 
Section 3.1). 

5. In Section 3.1, we obtained an estimate for the mean of A3, and 
we know that A4 n| + As n% = A3 n| . As before, we can com- 
pute the fractions 04 = % and as = % by linear regression 

"8 ,! 8 

(see Figure 4), and thus obtain an estimate for kg. Note that 
we will reject values of kg that imply the proliferation time is 
larger than 7 days (see Section 3.1). We note that the p-values 
for 03, fl4, and as are given by 8.43 x 10~ 3 , 3.33 x 10~ 7 , and 
4.56 x 10~ 8 , respectively (smaller than the significance level). 

6. From steps 3, 4, and 5 above, we have generated (a table of) val- 
ues for (pg,k4, and As, given a fiducial value for cp4 in the interval 



[0, <p2\. As discussed in Section 3.1, we only accept parameter 
values which provide thymocyte cell counts at time t = 21 days 
that are within one standard deviation from the experimentally 
determined values (see Table 1). 

We obtain the following parameter values: 

(14 = 0.04 day -1 , fi s = 0.11 day -1 , 
£4 = 0.21 day -1 , § 8 = 0.14 day -1 , 



and 



<P4 = 0.070 day 



: 0.054 day 



k 4 = 0.216 day -1 , k & = 0.093 day -1 



These parameters imply the following thymic selection rates: 

3.2.1. Death rates 

2.05 x 10 4 cells/h die by negative selection in compartment 4 
{1x4 n\) and 1.95 x 10 4 cells/h die by negative selection in 
compartment 8 (fig n% ). 

3.2.2. Differentiation rates 

2.50 x 10 4 cells/h are CD4 positively selected in compartment 2, 
that is, become CD4 SP from post-DP (<p 4 n|), 1.90 x 10 4 cells/h 
are CD8 positively selected in compartment 2, that is, become 
CD8 SP from post-DP {<p B n*), 1.08 x 10 5 cells/h leave compart- 
ment 4 to go to the periphery (£ 4 n|), and 2.48 x 10 4 cells/h leave 
compartment 8 to go to the periphery (fs «g). 



3.2.3. Proliferation rates 

11.10 x 10 4 cells/h proliferate in compartment 4 (A4 n\ 
1.60 x 10 4 cells/h proliferate in compartment 8 (kg «g). 



and 




10 20 30 40 



FIGURE 4 I Linear regression plots for the second model. 
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We can compute the stringency of thymic selection, defined by 
the ratio: 



g 4 n\ + & n\ 
0 



= 8.96% . 



We can provide an estimate for the cortical positive selection 
probabilities, that is the (per post-DP cell) probability to become a 
CD4 SP or a CD8 SP, and the probability to be negatively selected 
in the cortex. We have obtained: 



S4 



P2 = 



ip 4 



H2+<P4 + < 
M2 + <Pi + <Ps 



4.7%, s 8 



= 91.7% . 



M2 + <Pi + < 



3.6% , 



Finally, we have computed the (per cell) probability to die, given 
that the cell is in compartment as well as the (per cell) probability 
to proliferate in the medulla. We obtain: 



pi 



ix 4 



(I4 + £4 + A4 
/*8 + ?8 + ^8 



8.6% , <J 4 



= 32.1%, q s = 



k 4 



114 + £4 + X4 

A 8 

+ ?8 + ^8 



= 46.3% , 
= 27.0% . 



These probabilities imply that the probability to exit the thy- 
mus as a mature CD4 thymocyte (that has already reached the 
medulla) is given by 100 ~(^ 6 + 46 - 3 ) ; which is 45.1%, and the proba- 
bility to exit as a mature CD8 thymocyte (that has already reached 
the medulla) is given by 10Q ~ (3 1 2 0 q +27 -° ) , which is 40.9%. 

3.3. SENSITIVITY ANALYSIS 

In this section, we explore the sensitivity of the parameters to 
perturbations in the experimental data. For the first model, the 
experimental data are given in terms of the following eight 
quantities: 



0 = (Tu t 2 , T 3 ,0 out , fli, a 2 , «3, «i) 



where a\, a.2 are the coefficients of the linear regression of and 

ri* 

-|, respectively, and n, is the experimental mean value of n,-. 

"2 

We perturb each entry of the vector 9 by adding and sub- 
tracting 10% of its value. Therefore, we now have two values 
for 9j, equal to 9, + j^9j and 9, — jq9,. Consequently, we have 
2 8 sets of 9, which will be used to compute the corresponding 
model parameters as described in Section 3.1. Parameter val- 
ues will only be accepted if they provide a stable solution before 
t= 21 days. 

For the second model, the experimental data is given in terms 
of the following seven quantities: 

9 = (r 4 , r 8 , 1-L4, Ms. «3, fl4, <h) » 

where 03, a 4 , as are the coefficients of the linear regression of 
-I, -4, and -4, respectively. We have made use of the means of 
the following parameters of the first model: (p,(pi,fii,(p2,fJ-2>^3- 



Table 2 | Means, 95% trimmed and minimum-maximum intervals of 
the model parameters. 



Parameter 


Mean value 


95% Trimmed 


Minimum-maximum 






interval 


interval range 


<P 


35.86 x 10 6 


(35.65 x 10 6 , 


(28.93, 43.21 x 10 6 ) 




cells/day 


35.07 x 10 6 ) cells/day 


cells/day 


<Pt 


0.139 day" 1 


(0.138, 0.140) day" 1 


(0.112, 0.167) day" 1 


<P2 


0.136 day" 1 


(0.134, 0.139) day" 1 


(0.041, 0.274) day" 1 


<P4 


0.140 day" 1 


(0.136, 0.145) day" 1 


(0.060, 0.264) day" 1 


<P8 


0.134 day" 1 


(0.129, 0.138) day" 1 


(0.010, 0.214) day" 1 


M1 


0.265 day" 1 


(0.263, 0.267) day" 1 


(0.196, 0.333) day" 1 


V-2 


1.372 day" 1 


(1.365, 1.378) day" 1 


(1.083, 1.618) day" 1 


HA 


0.040 day -1 


n/a 


(0.036, 0.044) day" 1 




O.IIOday- 1 


n/a 


(0.099, 0.121) day" 1 


k 4 


0.181 day" 1 


(0.179, 0.184) day" 1 


(0.116, 0.226) day" 1 


>-8 


0.085 day" 1 


(0.080, 0.090) day" 1 


(0.078, 0.092) day" 1 


?4 


0.231 day" 1 


(0.230, 0.233) day" 1 


(0.229, 0.233) day" 1 


?8 


0.1 52 day" 1 


(0.150, 0.154) day" 1 


(0.149, 0.155) day" 1 



We perturb each entry of the vector 9 as described above. There- 
fore, we have 2 7 n, Pi sets of 9, with n m , the number of different 
values considered for q>4 in the interval (0, q>2). Parameter val- 
ues will only be accepted if they provide a stable solution before 
t= 21 days. 

The results of the sensitivity analysis, with 95% trimmed 
intervals 2 and minimum-maximum interval ranges, are given in 
Table 2. 

3.4. VARIABILITY IN THE SELECTION RATES 

The (trimmed and minimum-maximum) intervals derived in 
Section 3.3 allow us to estimate the variability in the different 
selection rates discussed in Sections 3.1 and 3.2. For example, 
given variations in the parameters, the corresponding variations 
in the selection rates can be shown to be: 



Api 
Ap 3 = 

A(J 3 = 

As,- 



1 



(Mi + <Pd 2 
1 



(M3 


+ <P3 + W 2 




1 


(M3 


+ <P3 + l-l) 2 




1 

r 



((PiAm + (ijA(pi) for i = 1,2, (10) 

[(</?3 + ^3) A/x 3 + M3 A^ 3 + M3 AA 3 ] , 

(11) 

[A 3 A^ + A 3 A<p 3 + (ju, 3 + <p 3 )AA 3 ] , 

(12) 

2 [<PiAfl 2 + <Pi&<Pj 



Api 



0*2 + <Pi + <Pj) 

+ O2 + <Pj)&<Pi] for i = 4, j = 8 or i = 8, j = 4 , 

(13) 

1 

: —j[fl,A$j + mAki 

+ (|j + Xi)AiXi\ for i = 4, 



(14) 



2 We define the 95% trimmed interval to be the result of the sensitivity analysis after 
trimming the lower and upper 2.5% of values. 
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(fii + & + A,) 2 

+ (& + /x,)AX,] for i = 4,8. (15) 
We present in Table 3 the variability of the selection rates. 
4. DISCUSSION 

We have brought together experimental data with a mathemati- 
cal compartment model [similar to other progression models of 
CD4 and CD 8 T cell development (13, 14, 18, 21, 33)] to provide 
estimates for the selection events that take place in the thymus. 
We have made use of a range of experimental data: (i) steady state 
thymocyte cell counts (26), mean residence times in each com- 
partment (27-29), murine thymic export rate (14, 26, 30), and 
recently reported asymmetric death rates for the CD4 SP and CD8 
SP thymocytes (21). Our preliminary results support the unex- 
pectedly high death rate in the post-DP thymocyte population 
observed in Ref. (21). We note that our approach is unrelated 
to that of Sinclair et al. both experimentally and mathematically 
(21). This rate, fi2> has been estimated to be at least an order of 
magnitude larger than any of the other death rates in the pre- 
DP, CD4 SP, or CD 8 SP pools (see Table 2). In terms of selection 
rates, our analysis yields the following: pre-selection thymocytes 
(pre-DPs) have a 65.8% probability of dying by neglect in the 
cortex, and a 34.2% probability of becoming post-selection thy- 
mocytes (post-DPs). At the post-selection stage, post-DPs have a 
91.7% probability of dying by negative selection (apoptosis) in 
the cortex, a 4.7% probability of becoming CD4 SPs, and a 3.6% 
probability of becoming CD8 SPs. In the medulla, CD4 SPs have 
an 8.6% probability of dying by negative selection (apoptosis), 
whereas CD8 SPs have a 32.1% probability of dying by negative 
selection. CD4 SPs have a 45.1% probability of exiting the thymus 
and reaching the periphery as mature thymocytes, whereas that 
probability for CD8 SPs is only 40.9%. Finally, the data supports 
some level of cellular proliferation in the medulla, with CD4 SPs 
having a 46.3% probability of proliferation and CD8 SPs a 27% 
probability. 

Earlier work by Mehr and collaborators combined experi- 
mental and theoretical approaches to estimate thymic selection 
rates (13, 33), neglected death rates in the medulla, but consid- 
ered potential feedback from mature T cells. In agreement with 
these authors, our results indicate that thymocyte death is high- 
est at the post-DP stage. However, as death in the medulla had 
been neglected, these authors concluded that the CD4:CD8 ratio 
in SP thymocytes is determined by the differentiation rates. In 
this paper, we have made use of CD4 and CD8, or medullary, 
death rates, which allowed us to directly compare cortical (DP) to 
medullary (SP) death rates. Furthermore, our approach allowed 
us to conclude that medullary, or SP, death was due to negative 
selection, as it was rescued by Bim deficiency (26). Sinclair et al. 
also recently addressed the temporal dynamics of thymic selection 
using an unrelated approach (both experimentally and mathe- 
matically) (21). While their experimental approach did not allow 
them to distinguish death by negative selection from death by 
other mechanisms, their overall finding was consistent with ours, 
that thymocyte death is highest at the post-DP stage. 



Table 3 | Selection rate values (initial and after perturbation) and their 
variability intervals. 



Rate 


Initial 
value (%) 


After 

perturbation {%) 


A Value 

(%) 


A Min— max 

(%) 


P1 


65.8 


65.66 


±0.76 


±20.55 


P2 


91.7 


90.98 


±0.49 


±17.28 


P3 


22.91 


24.07 


±0.90 


±28.22 


Q3 


42.24 


41.74 


±0.99 


±30.53 


S 4 


4.69 


8.51 


±0.80 


±15.16 


S8 


3.61 


8.13 


±0.79 


±15.03 


P4 


8.59 


8.84 


±0.37 


±4.91 


QA 


46.29 


40.07 


±1.29 


±20.49 


Ps 


32.12 


31.71 


±2.25 


±35.37 


Q8 


27.0 


24.5 


±3.58 


±64.81 



Further attempts to quantify thymic selection rates making use 
of mathematical models also include those of Faro et al. (12). The 
mathematical model developed by Faro and collaborators does 
not include time dynamics, but describes the relationship between 
the number of selecting ligands and the probability of selection 
of a given thymocyte. Thomas- Vaslin et al. (14) obtained esti- 
mates of thymic selection rates, using an experimental procedure 
that temporarily blocks thymic output and a mathematical model 
in which rates of transit from compartment to compartment 
depend on the number of cell divisions. Their model can capture 
the thymic "conveyor belt" (34, 35) scheme, but requires more 
differential equations and more parameters than equation (6). 
Despite the differences between their theoretical and experimen- 
tal models and ours, similar estimates for thymic selection rates 
are found. For example, we estimate that 1.2 million post-DP 
become CD4 SP thymocytes per day and 0.5 million post-DP 
become CD8 SP thymocytes per day; their estimates are 0.9 and 
0.2 million, respectively. Finally, we estimate that 2.6 million CD4 
SP thymocytes per day and 0.6 million CD8 SP thymocytes per 
day exit the thymus. Their estimates are 2.4 and 0.5 million, 
respectively. 

Our estimates of how many CD4 and CD8 SP thymocytes 
survive and exit the thymus reflect the skewed CD4:CD8 SP thy- 
mocyte ratio observed in C57BL/6 mice, which is approximately 
4:1 (36). This ratio is similar to the reported CD4:CD8 ratio of 
recent thymic emigrants (37), and raises the question of what 
accounts for the CD4 bias. While we were able to determine death 
and differentiation rates for both CD4 and CD8 SP thymocytes 
(see Table 2), our experimental approach did not allow us to deter- 
mine what fraction of the post-DP pool was MHC class I versus 
II restricted. Therefore, we could not address the issue of when 
and how the CD4:CD8 bias becomes established. The approach 
of Sinclair et al., which used MHC class I and class II deficiency, 
allowed them to address this question. Their data suggest that the 
skewed CD4:CD8 ratio reflects asymmetry in post-selection DP 
death rates, rather than more efficient positive selection of CD4 
compared to CD8 thymocytes (21). Yet, the parameter estimation 
allows us to compare the following different CD4:CD8 ratios (see 
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Section 3.2): (i) the CD4:CD8 ratio of positive selection in the 
post-DP pool (differentiation from post-DP to either CD4 SP or 
CD 8 SP) is given by ^ « 5 : 4, (ii) the CD4:CD8 ratio in the SP 

pool is given by -4 3 : 1, and (iii) the CD4:CD8 ratio of posi- 
tive selection in the SP pool (differentiation from SP to peripheral 
early thymic emigrants) is given by |p| ~ 4 : 1. Our observa- 
tions indicate that the CD4 bias is progressively established, as the 
thymocytes mature from the post-DP stage until the exit of the SP 
stage to migrate to the periphery. 

Our mathematical analysis has also allowed us to estimate the 
stringency of thymic selection, defined by: 

£4 til + £ 8 K 
= S8__8 = g 96% 

that is, the ratio between the number of thymocytes per unit time 
that exit the thymus and the number of thymocytes per unit time 
that enter the pre-DP stage. The sensitivity analysis described in 
Section 3.3 allows us to provide a value of Act = 0.2%, where 
we have made use of the minimum-maximum interval ranges 
(see fourth column of Table 2). A different measure of stringency 
could be based on the probability of a cell surviving the maturation 
process. In our notation, this would correspond to the following: 

(1 - pi) x (1 - p£) x (1 - = 2.19% . 



We note that this measure of stringency is the probability of 
not dying in any of the three compartments considered in the 
model (pre-DP, post-DP, and SP). As discussed in Appendix A.l, 
and given that in the SP pool, thymocytes may proliferate, there 
is a need to consider this special case. Our estimates suggest that 
a population of 10 3 pre-DP thymocytes will yield 69 CD4 and 25 
CD8 SP thymocytes that leave the medulla to get incorporated into 
the peripheral naive T cell pool (see details in Appendix A.l). 

The sensitivity analysis (see Section 3.3) and the variability of 
the selection rates derived from it (see Section 3.4) give us the con- 
fidence to conclude, that our parameter estimation is robust. We 
are aware that the experimental data we have made use of [steady 
state thymocyte cell counts (26) ] do not provide the exquisite time 
resolution described in Ref. (21). However, the supporting math- 
ematical model described in Section 2.2, allows us to obtain the 
time evolution of the thymocyte populations, once the parameters 
have been estimated. In Figure 5, we plot the time evolution of 
the total number of cells in each compartment of the mathemat- 
ical model: pre-DP, post-DP, CD4 SP, and CD8 SP thymocytes. 
We start with no cells at time zero, tii(t= 0) = 0 for i= 1, 2, 4, 8. 
Trajectories have been plotted for a period of 6 weeks and have 
been computed for every permutation of the parameter set pre- 
sented in Table 2. The subset of parameters shared with the simple 
model (</>, tpi, n\, /X2), were fixed at their mean values. Thus, 548 
distinct parameter sets were generated. The system of equations 
(6) was solved using a fourth order Runge-Kutta method (Python 
source code). 




FIGURE 5 |Time evolution of the thymocyte populations in the second model. The different trajectories correspond to the parameter values and ranges 

described in Table 2. 
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The approaches introduced in this paper have shed some 
light on the probabilities and timescales that characterize cellu- 
lar fate in the thymus after the DN stage. We plan to generalize 
the mathematical model introduced here, making use of experi- 
mental data for the strength of TCR binding in Nur77 GFP mice 
(26), to investigate issues such as the death rate in the post- 
DP pool and the CD4:CD8 ratio. Our model assumes that all 
progenitors in a particular pool behave with identical kinetics, 
i.e., move through the various stages of selection at the same 
rate. Future model refinements will come from consideration 
of the heterogeneity of the pools, which are known to include 
cells that will become iNKT cells, regulatory T cells, and intra- 
epithelial lymphocytes (2). It is also possible that progenitors of 
the same general class move through the selection process with 
different kinetics (34). The models introduced here can serve 
as a first step to study human thymic selection, although com- 
prehensive data on human thymic subsets, their sizes, and res- 
idence times are not yet available. It would be of great interest 
to apply the model to data on thymic subsets and cellularity in 
children, keeping in mind that residence times of human sub- 
sets may differ from murine ones (38). Finally, we note that we 
have not mentioned the relevance of cytokines, such as IL-7, 
during thymic development. Some differences have already been 
described for the role of IL-7R in human versus mouse T cell 
development (38, 39). We hope in the near future to combine 
mechanistic mathematical models of IL-7 and IL-7R (40) with 
the T cell development model introduced here to address these 
issues. 
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APPENDIX 

A.1. STRINGENCY OF THYMIC SELECTION: A STOCHASTIC MODEL 

In this section, we present the details that allow us to compute 
the stringency of thymic selection for the mathematical model 
considered in Section 2.2. 

Let us assume that at time t= 0, there exists a single T cell in 
a given compartment. In our case, the compartment can be the 
pre-DP, post-DP, or SP (either CD4 SP or CD 8 SP) stages. The 
cell, at any time, may die (with rate, /x), divide (with rate, X), and 
produce two daughter cells, or leave the compartment (with rate, 
£ ) to enter a different compartment. The waiting times for each 
event are assumed to be exponentially distributed, and daughter 
cells are assumed to behave identically to the initial single cell. 
We introduce the bivariate Markov process {X(t), Y(f)}f>o, where 
X(t) is the number of cells in the compartment at time f, and Y(t) 
is the number of cells which have left the compartment (to enter 
a different compartment). Our aim is to calculate the expected 
number of cells (and variance) that leave the compartment. 

State probabilities for the Markov process are defined as 
follows (41) 

p (x , y) (t) = Prob{X(0 = x, Y(t) = y\X(0) = 1, 7(0) = 0} , 



dmxx(t) 

dt 
dm X y(t) 

dt 
dmyy(t) 

dt 



2(X - ii - % )mxx(t) + (A - fJL - l)m x (t) , (A6) 
(X - /x - f)mx? 0) + %[mxx(t) - m x (t)] , (A7) 
f [2mxY(t) + m x (t)l (A8) 



Given that we start with a single cell, the expected number of 
cells at time t is given by 



m x (t) 



(A9) 



Under the restriction X < /x + £, the expected number of cells 
tends to zero as t— > +oo. This implies that all cells from the single 
T cell progenitor either die or leave the compartment for suffi- 
ciently large times. The expected number of cells which leave the 
compartment is given by 



my(t) 



+ £ - X 



h _ e 0.-M-f)t] 



(A10) 



(Al) ment can be shown to be 
and transition probabilities for this process are defined as fol- 
lows (41) 



As t — > +oo, the expected number of cells to leave the compart- 



lim my(f) 



(All) 



P(w,z),(x, 7 )(Af) = Prob{X(f + At) = w, 
Y(t + At) = z\X{t) = x, Y(t) = y} 

XxAt + o(Ar), if w = x + l,z = y , 

lix Af + o(At), if w = x — l,z = y , 

■ tjxAt + o(At), if w = x — l,z = /+ l, 

1 - (X + ii + %)xAt + o(At), if w = x,z = y , 
o(At), if w, z otherwise . 

(A2) 

The Kolmogorov (or master) equation for this process is given 
by (41) 



Ii + f - X 

We now solve the remaining ODEs equations (A6-A8), to find 



mYY(t) 



2X% 2 



{X — ii — £) 3 

4A£ 2 
~ (X - M - §) 2 



J.(l-ll-^)t _ (X-fi-^t 



1 



+ 



Jk-(i-i)t 



X- 11-% 

2X% 2 



JA.-^-f)f 



X- n-$ 
X- /u.-| 



{x-ii-^y 



(A12) 



dt 



X (x — l)p( x - + n (x + l)p (x+Uy) (t) 



+ f (X + l)p (x+ l,y-l)(t) - (X + II + f) X P(x,y)(t) . 

(A3) 

Let mx(t) be the expected number of cells in the compart- 
ment under consideration, and my{t) be the expected number 
of cells which have left the compartment. Similarly, let mxx(t) 
be the expectation of the random variable X(t) 2 , mxyit) be the 
expectation of the random variable X(t)Y(t), and myy(f) be the 
expectation of the random variable Y(t) 2 . Then, making use of 
the probability generating function technique (41), we derive the 
time evolution for the first two moments of the system: 



It, therefore, follows that the random variable Y(t), which 
represents the number of cells leaving the compartment under 
consideration, has the following variance (in the limit t — > +oo) 



<Ty = lim \niYY(t) 



my 



(t) 2 ] 



2X$ 2 



+ 



(fl + l- X) 



11 + %-X 
(A13) 



A.2. 



dm x {t) 

dt 
dmy(t) 

dt 



(X - n - l)m x (t) 



^mx(t) , 



(A4) 
(A5) 



STRINGENCY OF THYMIC SELECTION IN THE FOUR 
COMPARTMENT MODEL 

The previous example can easily (but laboriously) be extended to 
the mathematical model introduced in Section 2.2. 

We may evaluate the expected number of, for example, CD4+ 
T cells produced by a single pre-DP progenitor (or more gener- 
ally N pre-DP progenitors). We only present the time evolution 
of the moment generating function. The deterministic equations 
describing the mean and variance of the numbers of cells in each 
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compartment (pre-DP, post-DP, SP CD4, and SP CD8) are left for 
the reader to derive. When counting the number of CD4+ T cells 
leaving the thymus, the moment generating function satisfies the 
following partial differential equation 



3M n 3M a n 3M 



3M 

Of 2 



+ M4(e 



3M , 
1)— +A 4 (e^ 
9 6 4 



1) 



3M 

304 



1) 



3M 
304~ 



1) 



3M 

307 



(A14) 



The symmetry of the mathematical model implies that an 
equivalent equation for the number of CD8+ T cells leaving the 
thymus can be obtained by interchanging the indexes 4 and 8. 
For our derived parameter set, the previous equation allows us to 
conclude that the expected number of CD4+ T cells, a single thy- 
mocyte in the pre-DP compartment produces, is 0.069 (standard 
deviation 0.96), whereas the expected number of CD8 + T cells 
which leave the thymus is 0.025 (standard deviation 0.59). 

To put this into perspective, a population of 10 3 pre-DP thy- 
mocytes is expected to produce 69 CD4 + T cells which leave the 
thymus (standard deviation 66), and 25 CD8+ T cells (standard 
deviation 41). More generally, a population of N pre-DP thymo- 
cytes is expected to produce 0.069N CD4+ T cells and 0.025N 
CD8+ T cells. 
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